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ABSTRACT 

Modelling AGN feedback in numerical simulations is both technically and theoreti- 
cally challenging, with numerous approaches having been published in the literature. 
We present a study of five distinct approaches to modelling AGN feedback within 
gravitohydrodynamic simulations of major mergers of Milky Way-sized galaxies. To 
constrain differences to only be between AGN feedback models, all simulations start 
from the same initial conditions and use the same star formation algorithm. Most 
AGN feedback algorithms have five key aspects: the black hole accretion rate, energy 
feedback rate and method, particle accretion algorithm, black hole advection algo- 
rithm and black hole merger algorithm. All models follow different accretion histories, 
and in some cases, accretion rates differ by up to three orders of magnitude at any 
given time. We consider models with either thermal or kinetic feedback, with the 
associated energy deposited locally around the black hole. Each feedback algorithm 
modifies the region around the black hole to different extents, yielding gas densities 
and temperatures within r <~ 200 pc that differ by up to six orders of magnitude at 
any given time. The particle accretion algorithms usually maintain good agreement 
between the total mass accreted by Mdt and the total mass of gas particles removed 
from the simulation, although not all algorithms guarantee this to be true. The black 
hole advection algorithms dampen inappropriate dragging of the black holes by two- 
body interactions. Advecting the black hole a limited distance based upon local mass 
distributions has many desirably properties, such as avoiding large artificial jumps and 
allowing the possibility of the black hole remaining in a gas void. Lastly, two black 
holes instantly merge when given criteria are met, and we find a range of merger times 
for different criteria. This is important since the AGN feedback rate changes across 
the merger in a way that is dependent on the specific accretion algorithm used. Using 
the Mbh _ c relation as a diagnostic of the remnants yields three models that lie within 
the one-sigma scatter of the observed relation and two that fall below the expected 
relation. The wide variation in accretion behaviours of the models reinforces the fact 
that there remains much to be learnt about the evolution of galactic nuclei. 

Key words: black hole physics - galaxies: interactions - galaxies: active - methods: 
numerical 



1 INTRODUCTION 

In the hierarchical model of galaxy formation, the largest 
galaxies are formed last. Naively, we would expect the high- 
est star formation rates (SFRs) and the activity from Ac- 
tive Galactic Nuclei (AGNs) to occur in these most massive 
galaxies. However, observational evidence contradicts this, 
showing that in massive galaxies, the peak SFRs and peak 
AGN activity occurred at redshifts 1-2 (e.g. Shaver et al 



1996 Madau et al.p.9 96), and not today. This reduction in 
activity from z ~ 2 to today was termed 'downsizing' by 
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Cowie et al. (1996). One favoured explanation of downsiz- 



ing is that during mergers, gas from the merger fuels both 
star formation and AGN activity (e.g. Sanders et al.|[l988 



Scannapieco et al. 20051; the feedback from the increased 



AGN activity then blows away all the gas, leading to a red 
and dead galaxy (e.g. |Springel et al.|2005 |. 

The observational motivation for this picture is the ev- 
idence that a supermassive black hole resides at the cen- 
tre of all galaxies with stellar spheroids (e.g. |Kormendy| 



& Richstone 1995 Ferrarese & Merritt 20001, and that 



they did not evolve independently of one another. The two 
strongest correlations are the relationship between the black 
hole mass and the stellar velocity dispersion (Mbh-c; e.g. 
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Silk & Rees 1998, Ferrarese & Merritt 2000. Gcbhardt et al.| 


2000 


Tremaine et al. 2002 


King 2003 


|Giiltekin et al. 2009 


), 



and the black hole mass and the bulge mass (MeH-Mt,; 
Magorrian eTaT1[T998l |McLure fc Dunlop|[2002l |Marconi"fe"l 



Hunt 20031. A likely implication of these relations is that 
the AGN feedback is self-regulated: Outflows from the black 
hole following a strong accretion event interact with the sur- 
rounding gas, inhibiting further accretion events, and hence 
limiting black hole growth (e.g. [Silk 



Rees||1998| |Fabian 

1999b; Sc annapieco fc Oh|2004 l. Large X-ray cavities in the 



hot gas halo around AGN (e.g. Boehringer et al. 1993 Mc- 
Namara et al. 2000 1 suggest that some outflow events can 



be quite powerful; these cavities are likely formed by jets 
from the AGN (e.g. |McNamara fc Nulsen||2007[ ). Although 
there is significant evidence supporting the AGN scenario 
for explaining downsizing, it is not yet understood precisely 
how feedback energy couples to the surrounding gas; possi- 



ble mechanisms are line radiation pressure (e.g. Castor et al 



1975 


Proga et al. 


2000), radiation pressure on dust grains 


(e.g. 


Murray et al. 


2005 ) , Compton heating of Mailing gas 


(eg. 


Ciotti & Ostriker 


2001 1 and photo-ionisation pressure 


(Buff & McCray||1974 


Cowie et al.|1978l. 



AGN feedback has been implemented in many numeri- 



cal simulations (e.g. Kawata & Gibson 2005 Di Matteo et al. 


2005 


Springel et al. 2005 Thacker et al. 2006; Sijacki et al. 


2007 


Okamoto et al. 2008 


Kurosawa et al.||2009| |Booth & 


Schaye||2009| |Debuhr et al. 


2011 Dubois et al.||2012|), but 



both theoretical understanding and numerical implementa- 
tions are fraught with difficulties. First, the spatial scales 
involved in studying AGN and the related feedback span 
many orders of magnitude, and information on all of these 
scales is needed simultaneously to fully and properly under- 
stand AGN and their feedback. On the smallest scale is the 
Schwarzschild radius, 



2GMbh 

r 2 



(1) 



where Mbh is the mass of the black hole, G is Newton's 
gravitational constant and c is the speed of light; typical 
values range rs ~ 3 x 10 6 ~ 8 km. On the next larger scale is 



the Bondi radius ( Bondi|1952 |, 



fBondi = 



GMbh 



(2) 



where c x is the sound speed of gas at infinity. The value of 
the Bondi radius is dependent on the black hole's environ- 
ment and can range from a few parsecs to tens of parsecs 



(e.g. Kurosawa et al. 2009 Springel et al. 20051. This ra- 



dius, also known as the capture radius, divides a gas flow 



in to two regimes (Frank et al. 2002). Consider a spheri- 



cally symmetric gas cloud centred on a black hole, where 
the gas is initially at rest at infinity. The only forces acting 
on the gas are the gravitational force from the central black 
hole and the pressure forces within the gas (assuming we ne- 
glect the self-gravity of the gas). Beyond the Bondi radius, 
the gas is comparatively uninfluenced by the black hole and 
flows subsonically. Within raondi, the gas density begins to 
increase, and the gas flow inward eventually reaches a sonic 
point. At the sonic point, the gas plunges at a free-fall rate 
in to the black hole (iHobbs et al.||2012|). Thus, the Bondi 



radius can be viewed as the black hole's gravitational radius 
of influence. 



The last spatial scale of interest is that of an entire 
galaxy (or a galaxy cluster), which can span dozens of kilo- 
parsecs (or a few megaparsecs). When considering all of 
these scales, comparing the size of a black hole to its mas- 
sive host galaxy is similar to comparing a coin to the Earth 
( Fabian|20l"2| . 

To complement the range of spatial scales, studying 
AGN feedback also requires a large range of temporal scales. 
At short intervals, observations show the luminosity of the 
central engines of AGN varies on time-scales ranging from 



days to years (e.g. Webb & Malkan 2000 and references 
therein; [Sarajedini et al.||2011| ); moreover, there are short 
term differences in variability amongst the different classes 
of AGN, making the variable luminosity challenging to un- 
derstand and constrain. Large scale observations have de- 
tected large (~10 kpc radius) X-ray cavities in the gas 

around the AGN. To inflate these cavities, a n outburst of 
10 58_ 1Q 6i ergg of 

energy would be required (Birzan et al 



2004 



McNamara et al.|2005|) every few ~10' 

i 



averaged output of ~10 



10 45 ergs s" 



yr, or a time- 



I Churazov et al. 



20021. Thus, the next inherent difficulty in modelling AGN 
becomes obvious: AGN luminosity varies on time-scales as 
short as days, yet they are expected to produce major out- 
bursts every few ~10 8 yr. 

Numerically, we can draw a few parallels between AGN 
feedback and stellar feedback. Both have been added to nu- 
merical simulations to improve realism, and the result was 
better agreement with observations. However, star forma- 
tion and stellar feedback is a conceptually simpler problem 
than AGN feedback. In numerical simulations, if a packet of 
gas meets a given set of criteria, then stars are formed and 
feedback energy is returned; the star formation and feedback 
parameters can be reasonably well constrained using current 
observations, where there are detailed observations of out- 
flows at many different stages of stellar evolution, such as 
T-Tauri stars (e.g. |Cabrit et al.11990") jHartigan et al.|1995| , 
Wolf-Rayet stars (e.g. Crowther|2007 \ , mass loss and stellar 
winds from evolved stars (e.g. |Willson|2000| ) and supernovae 
(e.g. Reynolds 2008). AGN feedback is much harder to con- 
strain since there are fewer candidates for detailed observa- 
tion due to the peak of AGN activity being so far in the 
past. The nearest black hole candidate to study, Sagittarius 
A*, is currently in a quie scent phase; even durin g a recent 



(Nobukawa et al. 20111, its lumi- 



flare to 4 x 10 39 erg s" 1 
nosity remained belo w the typical AGN luminosity range 
of 10 40 -10 47 erg s" 1 (jFabianl 1999a|). In addition to uncer- 



tainty in modelling AGN feedback, tremendous amounts of 
energy can be numerically returned to the region around 
the back hole, creating steep gradients and presenting a 
challenge to numerical integration. Unlike numerical stellar 
feedback where there is difficulty in preventing unrealisti- 



cally rapid cooling after the feedback event (e.g. Gerritsen 
fc Icke||1997| |Mori et al.|[i"997| |Thacker fc Couchman||2000| " 



Springel & Hernquist 2003), continual AGN feedback may 
create shocks or unrealistically disrupt the system. 

When AGN feedback is included in numerical simu- 
lations, the approaches can vary widely, in both physics 
and numerical implementation. Also, the physics is often 
modelled using different numerical codes with different star 
formation algorithms and different initial conditions. Even 
starting from the same initial conditions, mass and spa- 
tial resolution, the different hydrodynamical approaches, gas 
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cooling, star formation and stellar feedback algorithms of 
different numerical codes produce different results, as shown 



in the comparison by Scannapieco et al. ( 2012 1; this compar- 



ison did not include AGN feedback. Currently, a direct nu- 
merical comparison of the AGN feedback mechanisms does 
not appear to have been published nor can it be compiled 
from the assembled literature. Here, we run four of the al- 
gorithms found in the literature using the numerical code 
Hydra. We also present a fifth simulation that is new to 
this study; this model is designed to take advantage of the 
'best' features of the other models. All of our simulations 
start from the same initial conditions and use the same star 
formation algorithm; this approach is conservative, but well 
constrained. We explicitly state that the goal of this is to 
highlight the different behaviours of the different algorithms 
and not to critique the various approaches. 

By using a pseudo-multiphase star formation model for 
all simulations, we unavoidably introduce some compromises 
in the AGN feedback models that rely upon a multiphase gas 
description; to compensate for this, we have introduced addi- 
tional variables to represent hot and cold fractions of the gas. 
Although this star formation model is not precisely equiv- 
alent to those implemented in other simulations of AGN 
feedback, we consider using it a necessary compromise to 
ensure that the variations in our results are only from the 
AGN feedback algorithm. 

The layout of this paper is as follows. In Section [2] we 
discuss our simulations, focusing on how the initial condi- 
tions are constructed and the AGN feedback algorithms. In 
Section[3] we will discuss our results, focusing on the impact 
of the different AGN feedback algorithms. In Section [4] we 
discuss the final state of each simulation, and we conclude 
with a review in Section [5] 



2 NUMERICAL SIMULATIONS 

To perform our simulations, we use the parallel version of 
Hydra ( Couchman et al.|1995 Thacker & Couchman 2006 1 , 
which uses an Adaptive Particle-Particle, Particle-Mesh al- 
gorithm ( Couchman||1991 1 to calculate gravitational forces 
and the standard smooth particle hydrodynamics method 
(SPH; |Gingold fc Monaghan|[l977l |Tucy|[T9T7| ) to calculate 
gas forces. It includes a star formation algorithm (see sec- 
and has been modified to include black holes and 



tion 



2.31 



AGN feedback (see section \2A\ 



2.1 Galaxy models 



To construct our model galaxy, we first use the GalactlCs 
package (Kuijken & Dubinski 1995 Widrow & Dubinski 



2005| |Widrow et al.| |2008[ > to create a Milky Way-sized 
galaxy that consists of a stellar bulge, stellar disc, and a 
dark matter halo; this is done through an iterative process 
to produce a self consistent system. The free parameters are 
chosen such that the component masses are similar to the 
component masses in Springel et al. (2005), and are given 
in Table [Q 



Component 


Parameter 


Value 


Bulge 




292 km s" 1 




Re 


0.7 lkpc 




n 


1.1 


Disc 


R i 


2.46 kpc 




~d 


0.49 kpc 




-Rtrunc 


30 kpc 




•Strunc 


1 kpc 




°"R0 


119 km s" 1 


Dark Matter Halo 


Oh 


13.6 kpc 




Hi 


275 kpc 




5r h 


25 kpc 






330 km s' 1 




7 


0.81 


metallicity 


Z 


0.05 Z 


mean molecular weight 


/' 


0.6 



Table 1. The chosen parameters for our model galaxies. All pa- 
rameters are defined in Section l2.1l of the text. 

The stellar bulge density profile is given by 



Pb(r) = p b 



Re 



-b(r/Ra) 1/n 



(3) 



which yields a Sersic law for the projected density pro- 
file if p = 1 - 0.6097/n + 0.05563/n 2 , where n is a 
free parameter. The constant pb is defined using Ob = 

|47rn6 n(p - 2) r [n(2 - p)] Rlp^^ , where al is the depth of 
the gravitational potential associated with the bulge and R c 
is the radial scale parameter, which is a free parameter; the 
variable b is adjusted such that R B encloses half the total 
projected light or mass. 

The stellar disc has a truncated density profile that falls 
off approximately exponentially in R and follows sech 2 in z\ 
the disc has radial and vertical scale heights Rd and Zd and 
truncation distances iitrunc and ztrunc- The radial velocity 
dispersion is given by a^(R) = <t^ &~ r ^ r " , where ctro is 
the central velocity dispersion and R a = Rd for simplicity. 

The dark matter halo profile follows 



Ph = 



V 2 



47ra 2 (r/a h )T(l + r/a h ) 3 - 7 



C(r;r h ,8r h ), (4) 



where ah is the halo scale length, rh is the cutoff radius, 
7 is the central cusp strength and o"h is a (line of sight) 
velocity scale that sets the mass of the halo. The truncation 
function, C(r; rh, 5rh) = |erfc ( J , smoothly goes from 
one to zero at r = rh over width Sr^- 

We then modify this initial galaxy in three ways. First, 
we convert ten per cent of the total stellar mass into gas to 
create a gas disc; the gas disc is the same as the stellar disc 
except that it has been reflected in the x — y plane to avoid 
coincidence with the star particles. Although the gas scale 
height is initially larger than physically motivated, cooling 
allows the gas to collapse into a thin disc within a few 10 
Myr; this produces a short transient evolution of the gas 
accompanied by a brief increase in the SFR. At resolutions 
higher than presented here, this vertical collapse produces a 
strong ring-shaped shock which propagates outwards. One 
solution is to relax the gas disc in a fixed potential and 
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then implement the relaxed disc in the initial conditions; 
another solution is to reduce the initial scale height of the 



gas disc. Studies in Williamson & Thacker (2012) show that 
at our fiducial resolution the pressure wave dissipates com- 
paratively quickly, and is largely insignificant at our lowest 
resolutions. Given also that we are interested in the rela- 
tive differences between the feedback algorithms during the 
merger, which is clearly a much more dynamic event, we 
consider the above gas disc construction acceptable. 

Second, we add a hot gas halo (hgh), which is chosen to 



follow the observationally motivated /3-profile (e.g. Cavaliere 
|fc Fusco-Femiano|1976| : 



Phgh(r) = po 



1 + 



(5) 



where po is the central density, r c is the core radius, and 
j3 is the outer slope parameter; we choose r c = 1.75 kpc 
and /3 — 2/3 as done in Moster et al. (20111. We set po by 



choosing the mass of the hot gas halo within 40 kpc to be 
equal to two per cent of the total disc mass ( jRasmussen et alT] 
2009); to conserve total halo mass generated by GalactlCs, 



we reduce the mass of the dark matter halo by the total mass 
of the hot gas halo. By assuming isotropy and hydrostatic 
equilibrium, the temperature profile of the hot gas halo is 
given by ( Kaufmann et aL||2007 l 

1 



kB Phgh(r) 



. .GMij). 
phgh(r) t. — dr, 



(6) 



where p, is the mean molecular weight, m p is the proton 
mass, kB is the Boltzmann constants, and M(r) is the to- 
tal mass interior to r. The hot gas halo is given an initial 
angular momentum which scales with the circular velocity, 
jhgh(ii) oc Rv c i IC (R), where R is the distance from the spin 
axis of the galaxy ( [Moster et al.|2011[ ). 

The final modification to the galaxy is the addition of 
a black hole (sink) particle, which is placed at the centre 
of mass and given an internal mass of 10 Mq. While this 
initial mass is more than ten times lower than expected from 
estimates of the mass of the Milk y Ways central b lack hole 
(M BH ~ (4.36 ± 0.42) x 10 6 M ; |Gillessen et al.||2009|>, it 
does match the initial black hole masses in Springcl et al.| 
( p005] > and |Debuhr et al.| poTTj ). Moreover, an initial mass 
lower than anticipated from the Mbh-o" relationship would 
be expected to grow quite quickly in a few Salpeter times 
because feedback at low masses is comparatively weak. At 
every step, we calculate a smoothing length for the black 
hole in order to calculate the gas properties around it, but 
the sink particle itself only experiences gravitational forces. 

Each galaxy has a total mass of 9.60 x 10 11 Mq, 1 287 
743 (168 351) particles for the fiducial (low) resolution sim- 
ulations, and a Plummer softening length of epiummcr = 120 
pc (epiummcr = 300 pc). The Plummer softening length 
is related to the S2 gravitational softening length, 632, by 
epiummcr = es2/2.34. See Table [2] for a breakdown of each 
galaxy. Lastly, we create a second, identical galaxy, and sep- 
arate them by 70 kpc; we then place them on parabolic tra- 
jectories around one another, in close agreement with the 



trajectories from Springel et al. (2005) 



2.2 Verification against other codes 

Since many of the models in the literature are run using 
Gadget2 ( |Springel et al.|2001| [Springel fc Hernquist|2002[ ), 
we ran the same adiabatic merger simulation with both Hy- 
dra and the publicly available version of Gadget2 to test 
the differences. Neither simulation included stellar or AGN 
feedback, thus we were essentially comparing the gravity 
and SPH solvers. The galaxies in each simulation followed 
the same trajectory, and synchronisation - as measured by 
the time to reach second periapsis - was within 0.2 per cent. 
Thus we can be confident that our results will be comparable 
with those currently found in the literature. 



2.3 Star formation 

The star formation algorithm implemented in Hydra is de- 



scribed and tested in Thacker & Couchman (20001; we will 



provide a brief summary here. This algorithm follows an 
approach to star formation that has now been studied ex- 
tensively in the literature, where star formation is allowed 
to proceed in regions where: 

(i) The gas exceeds the density limit of riH ~ 0.01 cm -3 , 

(ii) The flow is convergent, V ■ v < 0, 

(iii) The gas temperature is less than 3 x 10 4 K, 

(iv) The gas is partially self-gravitating: p g > 0.4pdm- 

When a gas particle meets the above criterion, the Schmidt 
Law (e.g. Katz 1992 Kennicutt||1998 1 is used to determine 
the amount of gas that is 'converted' in to stars. When the 
cumulative converted mass of a gas particle reaches half of 
the particle's original mass, m g , a star particle is spawned 
with mass m g /2 and the gas particle's mass is reduced to 
m g /2; when 80 per cent of the remaining gas mass is con- 
verted into stars, a second star particle is spawned with mass 
m g /2 and the gas particle is removed from the simulation. 
For computational efficiency, Hydra uses a Lagrangian form 
of the Schmidt Law, based upon the local density of the SPH 
particle, 



dM» 
At 



GrfrPg /2 M g 



where C slr is the star formation rate normalisation, p g is the 
SPH density of the gas particle, and M g is the mass of the 
gas in the particle that has not been converted in to stars. 
All other calculations in the code (e.g. density, smoothing 
length, etc..) assume that the entire particle is a gas particle 
of mass m g or m g /2. 

Whenever a star particle is created, feedback energy is 
immediately returned to the surrounding environment. Al- 
though our fiducial resolution models takes O(10 2 ) steps to 
evolve through the lifetime of an 8 Mq star, the lack of delay 
between star particle formation and feedback energy release 
is not a significant issue because at the resolution considered 
here we are still averaging over a number of giant molecu- 
lar clouds per particle and the evolution is representative 
at best. For every 100 Mq of stars formed, there is one 
supernova event, which contributes 10 ergs to the inter- 
stellar medium ( Sommer-Larsen et al.|1999 l, and feeds back 



5 x 10 15 e* ergs g of gas converted in stars, where e* is a 
dimensionless parameter; we set e* = 0.4 to match [Navarro| 



& White ( 1993 1 
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Fiducial resolution 



Low resolution 





Total 
(10 10 


mass 
M ) 


Particle mass 
(10 5 M ) 


Number of 
particles 


Particle mass 
(10 5 Mq) 


Number of 
particles 


Dark matter halo 




89.92 


11.75 


765 000 


89.92 


100 000 


Hot gas halo 




0.60 


0.36 


165 343 


2.77 


21 619 


Stellar bulge 




1.34 


2.37 


56 649 


18.10 


7 407 


Stellar disc 




3.56 


2.37 


150 375 


18.10 


19 662 


Gas disc 




0.54 


0.36 


150 375 


2.77 


19 662 


Black hole a 




10~ 5 


1.00 


1 


1.00 


1 



Table 2. Component breakdown for each galaxy. 

a Particlc (dynamical) mass is 10 9 M for Model DQM; internal mass of the black hole remains at 10 5 Mq. 



As first showed by Katz (19921, feedback energy re- 



turned to the interstellar medium (ISM) is radiated away 
quickly in high density regions. A reduced density is used in 
the cooling algorithm to prevent this immediate loss, thus 
allowing the feedback to influence the surrounding environ- 
ment; the reduced density decays back to its local SPH value 
with a half-life of tin = 5 Myr. The parameters e* and 



ti/2 are chosen to match those set in Thacker & Couchman 
( |2000[ > to reproduce the Milky Way's star formation rate in 
a simulation of an isolated Milky Way-like galaxy. 

This star formation algorithm is conceptually similar to 
the methods used in|Brook et al.| ( |2004| ), |Stinson et"aL] ( |2006| | 
and Governato et al. ( 2007 1, however, it varies from the sub- 



resolution multiphase models f ound in|Springel fc Hernquist| 
poM]), [Springel et al | ( [2005] ), |Booth fc Schaye| ( |2009| ) and 
Hopkins et al. (20091. In multiphase models, the subgrid 



physics is derived from a model of cold clouds (where stars 
form) embedded in a hot, pressure-confining phase. Above a 
given density threshold, the gas is thermally unstable to the 
onset of this two-phase medium. The mass fraction in each 
phase is determined by star formation and feedback, evap- 
oration of the cold clouds through thermal conduction, and 
the growth of clouds through radiative cooling. The star for- 
mation rate is then calculated based upon a prescribed law 
and an equation of state. A noteworthy difference between 
this method and our method is that here, gas can freely flow 
between phases, whereas in Hydra, once gas is 'converted' 
to stars, it is carried forward and not allowed to cool until 
the specified cooling period is reached. 



2.4 Black hole and AGN feedback algorithms 

AGN feedback algorithms essentially have five key compo- 
nents: 



(i) The accretion rate on to the black hole, 

(ii) The SPH particle accretion algorithm, 

(iii) The energy feedback algorithm, 

(iv) The black hole advection algorithm, and 

(v) The black hole merger algorithm. 

Each component will briefly be discussed below, then we will 
discuss the models in section [2~5l 



2.4-1 Accretion rates 

A commonly used accretion rate is the Bondi accretion rate 



(Bondi 19521 



Mr 



2nG 2 Mj HPoo 

(£4+^)3/2 



(8) 



where p x and Coo are the gas density and sound speed at 
infinity, v is the relative velocity between the gas at infinity 
and the black hole, and Mbh is the mass of the black hole. 
In hydrostatic equilibrium, the maximum physical accretion 
rate is the Eddington accretion rate, 

4nGM B nm p 



Me 



e r <r T c 



(9) 



where m p is the proton mass, <tt is the Thomson cross sec- 
tion, and e r is the radiative efficiency (i.e. the mass-to-energy 
conversion efficiency); we set e r = 0.1 (Sh akura &: Sunyaev 
|1973[ ). Typically, numerical simulations limit their accretion 
rate to the Eddington accretion rate, and we follow this con- 
vention. If the assumption of spherical symmetry is ignored, 
super-Eddington accretion rates can be achieved, which are 
associated with collimated outflows. In small scale outflow 



simulations by Kurosawa et al. ( 2009 1, they find steady-state 



results where super-Eddington accretion occurs in the equa- 
torial plane. 



2-4-2 Black hole mass growth and particle accretion 

In all cases, the 'internal' and 'dynamical' masses of the 
black hole is tracked. The internal mass, Mbh, is the mass 
of the black hole, which is increased by Msndt at every 
iteration; this mass is used in all calculations concerning 
AGN feedback. The dynamical mass, ttibh, is the mass of 
the sink particle, which is increased by the mass of a gas 
particle whenever one is accreted. The particle accretion al- 
gorithm should ideally maintain Mbh ~ tubh and directly 
address the loss of gas near the black hole due to accretion. 



In Model DQM (see [ 2.5.41, the dynamical mass is fixed for 
all time (except during a black hole merger), so the included 
particle accretion algorithm only simulates the loss of gas. 
When a particle is accreted on to the black hole, its mass 
and momentum are added to the black hole particle (except 
for Model DQM), and the gas particle is removed from all 
further computations; this accretion does not affect the in- 
ternal mass of the black hole. There are three categories of 
particle accretion algorithms: 
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(i) Stochastic- Unconditional: At all times, nearby parti- 
cles are tested to see if they will stochastically accrete. 

(ii) Stochastic-Conditional: Nearby particles are tested to 
see if they stochastically accrete if given criteria of the ttibh - 
Mbh relationship are satisfied. 

(iii) Continual- Conditional: Nearby particles are continu- 
ally accreted while given criteria of the ?tibh-Mbh relation- 
ship are satisfied. 



2.4.3 Feedback 

In all simulation, a portion of the accreted mass is returned 
as feedback energy, E = eMBHC 2 dt, where e is a dimen- 
sionless efficiency parameter. This energy is returned to the 



2.5 The Models 



nearby particles (except Model ONB; see { 2.5.3 1 either by 
increasing their internal energy (thermal feedback models), 
or by increasing their momentum via p = E/c (kinetic feed- 
back models). 



2.4-4 Black hole advection 

Properly tracking the black hole particle is critical since ac- 
cretion rates (hence feedback rates) depend on the local gas 
properties around the black hole. When the black hole mass 
is similar to a star or gas particle mass, it can be easily and 
inappropriately dragged around by two-body forces lead- 
ing to the inaccurate calculation of gas properties. Avoiding 
such behaviour is clearly desirable. In all of the simulations, 
a black hole advection algorithm is implemented to min- 
imise the inappropriate movement; except in Model DQM 



(see [ 2.5.4 1, the black hole's position is artificially updated 
after the completion of the gravitational solver algorithm 
but before the calculation of the accretion rate. 



2.4.5 Black hole mergers 

The process by which black holes merge is still a matter of 



active research (e.g. Escala et al. 2004 Berentzen et al. 2009 



Khan et al.|2011| |Bode et al.||2012[ ). Given that our simula- 
tions result in a combined halo and stellar system producing 
drag, it is reasonable to expect that the black holes would 
merge. Since any merger occurs on sub-resolution scales, 
models include a merger prescription to instantly merge the 
black holes when given criteria are met. When two black 
holes merge, one sink particle is removed and the remaining 
sink particle has the combined mass (both internal and dy- 
namical) and momentum, and is repositioned to the centre 
of mass of the two progenitors. 



2.4-6 Black hole's local environment 

Around every black hole we define a radius of influence, 
rinf. All gas particles within n n f contribute to the accretion 
properties at the black hole, are eligible to receive feedback 
energy (except for Model ONB), and are eligible to be ac- 
creted on to the black hole particle. Gas particles outside of 
ri n f have no explicit impact on the black hole or the AGN 
feedback algorithms. We set rj n f = max(2ftBH, 2ft m i n ), where 
a sphere with radius 2/ibh around the black hole particle in- 
cludes 60 gas particles, and h m i n is the smallest resolved 
smoothing length in the SPH solver. 



In sections |2.5.1| to |2.5.5| we describe the particular AGN 
feedback algorithms of our five primary models; a summary 
of each model can be found in Table [3] We ran four addi- 
tional models, each of which is a slight variant of a primary 
model; these are described in section [2.5.6| Every model was 
run at both fiducial and low resolutions. For nomenclature, 
we name these models after the initials of the authors of the 
paper it originally appeared in. The variant models have the 
same name as their parent model, followed by a lower case 
character to signify the difference. A model name followed 
by a subscript '1' explicitly refers to the low resolution ver- 
sion of the model; no subscript will refer to either the fiducial 
resolution version or to both versions, depending on context. 



2.5.1 Model SDH 



This model is based upon that found in Springel et al. ( 2005 1 



(herein SDH05). The accretion rate is given by a modified 
Bondi accretion rate, 



M B = 



4iraG 2 M^hP 



3/2 • 



(10) 



where c s and p are the local sound speed and density of 
the gas, and v le \ is the relative velocity of the black hole to 
the nearby gas. The free parameter, a, is included to relate 
the numerically calculated gas density and sound speed to 



what one would expect in reality. Booth & Schaye ( 2009 1 



argue that modest resolutions underestimate the density and 
overestimate the sound speed by orders of magnitude, thus 
justifying large values of a. As in SDH05, we set a = 100. 
Finally, the accretion rate is Eddington-limited, thus Mbh = 

min (Mb, M Ed d) . 

A given fraction of the accreted mass is allowed to re- 
turn to the surrounding environment as feedback energy. 
The rate of return is 



Ef c 



eie r M B nc 



(11) 



where e r = 0.1 is the radiative efficiency, and ef — 0.05 is the 
fraction of energy that can couple with the gas. This energy 
is returned kernel- weighted to the gas particles within n n {. 

To track the growth of the black hole particles, 
a stochastic-unconditional particle accretion algorithm is 
used. At every iteration, a probability, pi, is calculated for 
each particle, i, within r^f. This probability is then com- 
pared to a random number, Xi, and if pi > Xi, then the par- 
ticle is accreted. In this algorithm, the probability is given 
by 



Pi = WiMsup dt, 



(12) 



where Wi is the kernel weight of gas particle i relative to 
the black hole, and Xi is uniformly distributed on the the 
interval (0,1). 

To minimize inappropriate motions of the black hole 
particles, at every iteration the black hole is relocated to 
the nearby gas particle with the minimum potential energy 
provided that v IC \ < 0.25c s . If no gas particle meeting the 
velocity criteria exists within r lrl i, then the black hole is not 
advected. Once ttibh > 10m g , the black hole advection is 
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Model SDH 


Model BS 


Model ONB 


Model DQM 


Model WT 


Accretion Rate 


Mr (« = 100) 


M B (a ee a(n H )) 


Mdrag 


M visc 


M B (a = 100) 


Energy Feedback 


E = O.Q05M BH c 2 


E = Q.015M BH c 2 


E = 0.1L jct 


P = IQL/c 


E = 0.005M BH c 2 


FB Distribution 


kernel- weighted 


random packets of E CT ^ 


to 40 low-p particles 


isotropic 




Particle Accretion 


stochastic-unconditiona 


stochastic- conditional 


continual-conditional 


continual-conditional 


continual-conditional 


BH Advcction 


gas particle with 


gas particle with 


Ai ONB 




AZ WT 




v rcl < 0.25c B 


with v rel < 0.25c 8 


along stellar gradients 




towards centre of mass 




and smallest U 


and smallest U 








BH Merger 


d < h BH 


d < h BH 


d < ecj2 


d < h BH 


d < h BH 






& ^rcl < ^circ 


& gravitationally bound 




& "rel < c s 



Table 3. A summary of our five primary models describing each key algorithm. More detail and definitions of the variables are given in 
subsections |2, 5. l] to |2. 5.5] Numerically, /ibh = niax(/iBHi ^min)> where a sphere with radius 2/ibh around the black hole particle includes 
60 gas particles, and /i m i n is the smallest resolved smoothing length in the SPH solver. 



turned off and its movement is handled only by the gravita- 
tional solver. 

Lastly, two black holes merge when they come within 
each other's smoothing lengths and have a relative velocity 
less than the local sound speed; SDH05 argue that the local 
sound speed represents a simple measure of the velocity scale 
at which the black holes are able to merge. 



2.5.2 Model BS 



This model is based upon that found in |Booth fc Schaye] 
(20091 (herein BS09); this model has many similarities to 



Model SDH, but was originally implemented in a cosmolog- 
ical simulation. We caution that implementing their model 
in a higher-resolution simulation could lead to unwanted be- 
haviours as a number of the model parameters were chosen 
for their specific resolution. This model uses the (modified) 



Bondi accretion rate given in ( 10 1. If the resolution is suffi- 



cient enough, BS09 argue that the justifications for a = 100 
given in section |2.5.1| break down in low-density regions. 
Thus, this model sets a to be a function of the local hydro- 
gen density, nn- 



1 if riH < jih 
( "f I otherwise 



(13) 



where is the critical value required for the formation 
of a cold interstellar gas phase, and /3 is a free parame- 
ter; as in BS09, we set = 0.1 cm -3 and ,5 = 1. Fi- 
nally, the accretion rate is Eddington-limited, thus Mbh = 
min(MB,M Edd ). 



The rate of feedback is also given by (111, with e r = 0.1 
and ef = 0.15, which is three times more efficient than Model 
SDH. The feedback energy is allowed to accumulate until 
E > Edit, at which point a random gas particle within ri n f 
receives E CI it energy; this is repeated until the accumulated 
energy drops below -E cr i t . The critical energy is defined as 

m g fc B Ar 



-Ecr 



(14) 



(7 - 1) /wtih ' 

where m g is the (initial) mass of a gas particle and AT 
is the temperature increase a particle experiences with ev- 
ery feedback event. Due to the higher resolution in our 
model (m gas , fiducial = 3.6 x 10 4 Mq compared to 8.64 x 10 r 
Mq in BS09; e n duciai = 120 pc compared to 2 kpc h~ x 
after 2 = 2.91 in BS09), we set AT = 5 x 10 6 K as opposed 



to AT = 10 8 K used in BS09. Using the value in BS09 pro- 
duces very large temperature gradients in the SPH solver 
that are difficult to integrate accurately. 

Gas particles are accreted by a stochastic-conditional 
particle accretion algorithm. If Mbh < tibh, then the prob- 
ability of accretion is p, = 0, otherwise it is calculated using 



Pi = Wi (Mbh — wbe) p 



(15) 



As in Model SDH, particle i is accreted if pi > x%. 

The black hole advection is the same as in Model SDH. 
Lastly, two black holes merge when they come within each 
other's smoothing lengths and have a relative velocity less 
than the circular velocity at the radius of the most massive 
black hole's smoothing length. BS09 explicitly state that 
this merging criteria differs from SDH05 since the feedback 
returned may temporarily increase the local sound speed, 
thus may not be a representative velocity scale. 



2.5.3 Model ONB 



This model is based upon the model found in Okamoto ct al. 
( |2008[ > (herein ONB08), which was originally implemented 
in a simulation with cosmological initial conditions. We also 
note that this model is distinctly different from our other 
models in that it is specifically designed to reproduce the ra- 
dio mode of feedback. In this model, it is assumed that radi- 
ation from stars (e.g. through starbursts or winds) interacts 
with the rotating, clumpy ISM. This radiation irradiates one 
layer of gas at a time, extracting angular momentum from 
it. This permits an inflow of gas towards the galactic centre, 
and ultimately on to the black hole itself (jUmemura et al.| 



1997| |Umemura 2001 K awakatu fc Umemura||2002[ ). ThusT 
an accretion rate of gas on to the black hole can be calcu- 
lated by considering the stellar clouds in the region of star 
formation (RSF) near the black hole. Using these assump- 
tions, ONB08 calculate the drag accretion to be 

Lb,S¥ 



Ma 



cdrag 



')■ 



(16) 



where £dr 



1 is the drag efficiency, Lrsf is the total 



bolometric luminosity of all the stars in the RSF, and trsf 
is the total optical depth of the RSF. The total luminos- 
ity is calculated by summing the age-dependent bolomet- 
ric luminosities, which are obtained from a lookup table 
generated by pegase2 ( |Fioc fc Rocca-Volmerange||1997[ ). 
Next, the total optical depth of a gas cloud is given by 
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T c = XdPcTc — Xdm c /r'i where p c , m c and r c are the density, 
mass and radius of the cloud, respectively. Assuming that 
all the clouds are identical and randomly distributed over 
the region of star formation, the total optical depth can be 
approximated by 

3xd M c 



TRSF = 



(17) 



where \d = 50 cm 2 g _1 is the mass extinction coefficient, 
M c is the total mass of the clouds in the RSF, and -Rrsf is 
the radius of the RSF; since Hydra does not explicitly track 
a multiphase gas, we set M c to be half of the total gas mass 
within the RSF. We initialise -Rrsf = max (-R40, 2/i m i n ), 
where a sphere with radius -R40 centred on the black hole 
contains 40 gas particles. Then, larger radii are searched to 
find the radius that maximises Md ra g. Lastly, if the gas den- 
sity within the sphere is less than pthrosh = 5 x 10 -25 g cm~ 3 , 
the accretion rate is set to zero. 

In this model, it is explicitly assumed that the feed- 
back heats the halo gas through the production of jets. 
The jet mechanism used here is that of Meier (20011, and 



generates power from the rotational energy of the accre- 
tion flow and from the black hole itself. The accretion flow 
is divided into two regimes: standard thin discs (SD; op- 
tically thick, geometrically thin, radiatively efficient) and 
radiatively inefficient accretion flows (RAIF; optically thin, 
geometrically thick). Using the parameters in ONB08, the 
respective accretion-dependent luminosities are 

r SD 



r RAIF 

-'jet 



8.1 x 10~ 5 M BH c 2 
2.6 x IO^Mrhc 2 



if m > m cr it (18) 
if m < m crit , (19) 



where m = Men/A-fedd and m cr it ~ a 2 is the critical accre- 
tion rate that sets the division between the SD and RAIF 
regimes (e., 



Narayan et al. 19981. As in ONB08, we set 



Qsd = oraif = a = 0.1, where «sd an< A oraif are the 
viscosity parameters for SD and RAIF's, respectively. The 
feedback rate is then given by -Efeed = £r£jet, and the energy 
is distributed equally to the 40 nearest diffuse gas particles 

With p < O.lpthrcsh- 

To accrete particles, this model uses a continual- 
conditional particle accretion algorithm; when the internal 
black hole mass exceeds its dynamical mass, nearby gas par- 
ticles are accreted with the probability 

M B h - m B H 

P* = AT ' 



(20) 



where -/Vrsf is the number of gas particles within the region 
of star formation, and rrii is the mass of the i'th gas particle; 
particles are accreted until the dynamical mass exceeds the 
internal mass. 

To track the black hole, at every iteration, the local stel- 
lar density fields are computed and the black hole is moved 
along the steepest gradient by an amount 



AIqnb = min(0.01e S 2, 0.03 \v\ At), 



(21) 



where es2 is the gravitational softening length, \v\ is the 
velocity of the black hole, and At is the time-step; these 
coefficients are the same as in ONB08 and were determined 
empirically. 

Lastly, two black holes are assumed to merge when they 



are within a softening length of one another and are gravi- 
tationally bound. 



2.5.4 Model DQM 



This model is based upon that found in Dcbuhr et al. (2011 1 



(herein DQM11), and models accretion through the trans- 
portation of angular momentum, which is based upon multi- 
scale SPH simulations by Hopkins & Quataert (20101. The 
accretion rate is 



M v 



3tt<5E 



(22) 



where 5 is the dimensionless viscosity, E is the mean gas 
surface density, and f2 = \J GA'I/rf n{ is the rotational angular 
velocity of the gas. DQM11 treat 8 as a free parameter, and 
it characterises both the efficiency of angular momentum 
transport and the fraction of gas that is being converted 
into stars versus being accreted on to the black hole; as in 
DQM11, we set S = 0.05. 

The feedback is returned as momentum using 



L 

c 



(23) 



where r is the infrared optical depth, and L — 
min ^e r M v i ac c 2 , 1/Edd j is the Eddington-limited luminosity. 
This momentum outflow is used to approximate radia- 
tion pressure produced by absorption and scattering of the 
AGN's feedback; specifically, the ultraviolet radiation (emit- 
ted from the black hole) will deposit puv = L/c on to the 
gas, while infrared radiation (re-emitted from dust) will de- 
posit pir = tL/c on to the gas. Thus, the total p can be 
approximated as (1 + t)L/c~ tL/c for r > 1, which is valid 
near the peak of AGN activity when the black hole gains 
most of its mass. As in DQM11, we set r = 10. Lastly, the 
momentum is returned radially, such that every gas particle 
within ri n f receives an equal acceleration. 

In this model, there is no explicit artificial black hole ad- 
vection algorithm. Instead, a tracer mass is used to represent 
the black hole particle; that is, the dynamical mass of the 
black hole is initialised to »tibh = 10 9 Mq and held constant 
throughout the simulation. The black hole particle is now 
only advected by the gravitational solver; since the tracer 
mass in the fiducial resolution simulations is ~10 3 m,DM and 
~3 x 10 4 m g , it will not undergo artificial dragging by the 
surrounding particles. By initialising ztzbh = 10 9 Mq rather 
than 771bh = Mbh, the mass of each galaxy is increased by 
0.01 per cent. This means that Model DQM has slightly (but 
unavoidably) different initial conditions than the remainder 
of the models. The internal mass of the black hole is still 
initialised to 10 s Mq , and this mass allowed to grow as cal- 



culated by the accretion rate given in (221 



Since the dynamical mass of the black hole is fixed, a 
continual-conditional particle accretion algorithm removes 
'accreted' gas particles, but does not add their properties to 
the dynamical mass. Here, a random particle within rinf is 
removed whenever there is a mismatch between the amount 
of gas accreted via ( |22| and the total mass of removed gas 
particles. 

Lastly, black holes merge when they come within ri n f of 
each other, regardless of velocity. In this model, black hole 
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mergers are the only mechanism to increase the dynamical 
mass of the black hole. 



2.5.5 Model WT: This study 

This model uses the modified Bondi accretion rate given in 
( 10 1 for both its simplicity and for its wide use in the liter- 
ature (e.g. SDH05; |Robertson et al.|2006| [Croft et al.|2009 



Johansson et al. 20091; as in Model SDH, we set a — 100. 
The energy feedback rate is given by (111, and the feed- 



back energy is equally distributed amongst all the gas par- 
ticles within ri n f. This provides an isotropic heating to the 
core rather than preferentially heating (or super-heating) 
the particles very near the black hole, which may be tran- 
sient. 

The black hole advection algorithm is a modified ver- 
sion of that presented in Model ONB. First, the black hole 
is displaced towards the centre of mass of the sphere with 
radius r^nf centred on the black hole rather than along stellar 
gradients. This method still gives preference to the stellar 
distribution, but also considers the dark matter and gas dis- 
tributions. Second, the distance the black hole is displaced 
has been modified to 



AZwt = min(0.10/i B H, 0.30 \v\ dt, d CM ) 



(24) 



where dcM is the distance from the black hole to the centre 
of mass. We choose to use the black hole smoothing length 
rather than the softening length since all of the properties 
near the black hole are calculated using 2/ibh as the char- 



acteristic length. As in (21 1, our coefficients are empirically 



chosen so that non-negligible displacement would be possi- 
ble. Specifically, for the first option, /ibh < £S2, so its coeffi- 
cient needs to be increased; for the second option, our reso- 
lution is higher than in ONB08, thus we will have a smaller 
dt, thus we again need a larger coefficient. We choose this 
method rather than potential-well method of Models SDH 
and BS for the following reasons: If the AGN feedback cre- 
ates a gas void around the black hole, then either there will 
be no artificial displacement (allowing for the possibility of 
unnatural motion) or the black hole will be coupled to a gas 
particle on the edge of the void. In this method, as in Model 
ONB, the black hole only moves towards a particle rather 
than being coupled to it. 

As in Models ONB and DQM, we include a continual- 
conditional particle accretion algorithm: When A/bh > 
7?ibh + m g /2, we accrete the gas particle that is nearest 
to the black hole. The term m g /2 forces the dynamical and 
internal mass to oscillate around one another, and choos- 
ing the nearest particle is to reduce random effects in the 
simulation. 

The black hole merger algorithm is the same as given 
in Model SDH. We choose this algorithm since it had a dis- 
tance and a velocity criteria, and its velocity requirement 
is less stringent than in Models BS and ONB. The merger 
algorithms in Models BS and ONB were created for cosmo- 
logical models, whose resolution is lower than our fiducial 
resolution runs, thus we argue that our implementation of 
their merger prescriptions are more stringent than they in- 
tended. 



2.5.6 Additional models 

We have tested four additional models, each of which is a 
slight variation of one the above models. These models are 
identical to their parent models described above, except for 
the variation listed below. 

• Model BSw: This uses the black hole advection algo- 
rithm of Model WT; early tests of Model BS showed very 
erratic black hole motion that could possibly compromise 
the results. 

• Model ON Be: This uses a very conservative search al- 
gorithm to calculate Mdra g ; this yields a very small -Rrsp 
and a very low accretion rate. 

• Model WTh: For a resolution test, this model uses 
ft-min — > /imin/2; this should only impact calculations per- 
formed in very dense regions. 

• Model DQMe: This model uses r- m f = 4es2 = 1.17 kpc 
(compared to r inf ~ 73 pc of Model DQM); in DQM11, 
they fix rinf = 4es2, although this value is ~188 pc in their 
models. 



3 RESULTS 

Each of our models was evolved through a merger event, 
similar to that of SDH05. Model ONB was evolved for 1.1 
Gyr and the rest were evolved for 1.5 Gyr (including Model 
ONBi). By returning the feedback energy to the halo gas 
in Model ONB, a dense galactic core formed near the core 
merger epoch, which resulted in an extremely large wall- 
clock time per step for the fiducial resolution version. 

Each model followed a similar qualitative history, which 
is shown for Model WT in Figs. [I] and [2] for the gas column 
density and gas temperature, respectively. From this evolu- 
tion, we see four significant epochs: first periapsis, apoapsis, 
second periapsis, and core merger. The times of these epochs 
are calculated using the black holes as proxies for the centres 
of the galaxies, thus the epochs essentially represent the local 
minimum and maximum separation of the two black holes. 
In all models, first periapsis occurs at 166 Myr and apoapsis 
occurs at 480 Myr; the latter is sustained for approximately 
100 Myr. Second periapsis occurs approximately 11 Myr ear- 
lier for Model DQM than for the rest of the models. We have 
verified that this more rapid evolution is due to the inclusion 
of a more massive sink particle representing the black hole. 
In two-body simulations, where each particle represents one 
galaxy and dynamical friction is avoided, the period of the 
particles representing galaxies in Model DQM is 3.4 Myr 
(0.23 per cent) shorter than the particles representing the 
galaxies without tracer masses. For the full system, we find 
that the maximum separation at apoapsis is very slightly 
smaller (1.1 per cent), likely due to slightly higher dynami- 
cal friction, leading to an earlier second periapsis. Lastly, the 
core merger takes approximately 10 Myr to complete, and 
the onset of this process spans a range of 15 Myr amongst 
the different models. See Table [4] for a list of when these 
epochs occur. 



As might be expected, there are notable morphologi- 
cal differences between the models; this is readily apparent 
in Fig. [3] which displays the gas temperature, gas column 
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Figure 1. Evolution of gas column density for Model WT. Times from the onset of the simulation are listed in each frame; each frame 
measures 100 kpc per side, with an image resolution of 98 pc pixel . 





Dynamic Mass 


Tracer Mass 


First Periapsis/Gyr 
Apoapsis/Gyr 
Second Periapsis/Gyr 
Core Merger/Gyr 


0.166 (0.166) 
0.480 (0.480) 
0.884 (0.870) 
0.987 (0.973) 


0.166 (0.166) 
0.480 (0.480) 
0.872 (0.861) 
0.970 (0.962) 



Table 4. Important epochs from the onset of the simulation for 
the fiducial (low) resolution models. Average times are listed when 
a range is available. The Dynamic Mass models are Models SDH, 
BS, ONB and WT; the Tracer Mass model is Model DQM. 



density and stellar column density of the top left galaxy in 
each model at apoapsis. By apoapsis, a bar has developed 
in our non-tracer mass models (Models SDH, BS, ONB and 
WT), with bar strength depending on model. We have veri- 
fied that the lack of bar formation in Model DQM is a result 



of the additional mass from the tracer mass and not from the 
kinetic feedback. This extra mass causes an increase in the 
rotational velocity curve of the galaxy. See Figure [4] where 
we have plotted the (S2 softened) rotation curve for five 
variations on Model WTi just prior to the onset of bar for- 
mation; in this test suite, a tracer mass represents the black 
hole and has been varied between 10 6 and 10 9 Men. The ro- 



©■ 

= 10 e 8 Mq are similar, and the 



tation curves for M tl 
peak velocity is 2.5 per cent higher for M tr accr = 5 x 10 B 
Mq and 4.0 per cent higher for M traC cr = 10 9 Mq. The 
higher rotational velocities for M traccr = 10 9 Mq stabalises 
the galaxy against bar formation. 



To ensure that our initial black hole mass is reasonable, 
we ran Model WTi with seed black hole masses of 10 5 , 10 6 
and 10 7 Mq. All three models had final black hole masses 
within 1.5 per cent of one another, and all three remnants 
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Figure 2. Evolution of gas temperature for Model WT. Times from the onset of the simulation are listed in each frame; each frame 
measures 100 kpc per side, with an image resolution of 98 pc pixel -1 . The top right panel at 0.23 Gyr has a large filling factor of 
supernova events, caused by tidal interactions at first periapsis. Although ubiquitously occurring through the galaxies, the global star 
formation rate is at a local minima (see section [3. 6|. 



had A4bh~o" relations that agreed with the observed relation 
within the one-sigma standard deviation. 



3.1 Black hole advection 

The location of the black hole with respect to the galactic 
core plays a fundamental roll in determining the gas prop- 
erties used to calculate the accretion and energy feedback 
rates. Even small displacements with respect to the centre of 
the potential well can have notable effects. We have studied 
four artificial advection algorithms, each yielding different 
results. Since the results presented in the rest of this pa- 
per are implicitly coupled to the behaviour of the advection 
algorithm, we discuss this issue first. See Fig. [5] where we 



have isolated the path of one black hole in Models SDH, 
ONB and DQM. 

Without any artificial advection, low mass black hole 
particles are easily dragged by the local particles (whose 
mass may be comparable or greater than the mass of the 
black hole); even with artificial advection, some black hole 
'chaotic' motions (or oscillations) persist, depending on the 
algorithm. The smoothest path is from the tracer mass 
(Model DQM); however, the increased mass of the galaxy 
slightly adjusts its trajectory, decreases its evolution time 
compared to the rest of the models, and impacts the disc 
morphology (i.e. prevents bar formation). The least smooth 
path is from coupling the black hole to the gas particle with 
the minimum potential energy, provided that v IC \ < 0.25c s 
(Models SDH and BS). Frequently, there is no gas particle 
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Figure 3. Top loft galaxy of each model at apoapsis. Varying morphologies for the different models emerge as the galactic discs reform 
after first periapsis. The lack of bar in Model DQM is a result of the more massive black hole particle (i.e. the tracer mass). Each frame 
measures 20 kpc per side, with an image resolution of 39 pc pixel . From left to right: Gas temperature, gas column density and stellar 
column density. 



A Comparative Study of AGN Feedback Algorithms 13 



250 




0123456789 10 
Radius (kpc) 

Figure 4. The rotation curves for five low resolution test mod- 
els, where the mass of the tracer mass has been varied. For each 
model, the rotation curve is for the top left galaxy, and is calcu- 
lated at 0.25 Gyr just prior to the onset of bar formation. The 
parent model is Model WTj. 
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Figure 5. A portion of the unaveraged black hole path (including 
the turnaround point at apoapsis) in Models SDH, ONB and 
DQM for the fiducial (left) and low resolution (right) simulations. 
The dot indicates where the black hole advection algorithm turns 
off for Model SDH. The origin of the system is (x, y, z) = (0, 0, 0) 
and the black hole originates at (—27.6, —21.9,0.0). 



satisfying the velocity criterion, thus the black hole oscillates 
about the galactic core. When there is a gas particle satis- 
fying the criteria, it is likely near the periphery of n n f since 
the inner particles are typically too hot. Thus, the stringent 
velocity criterion permits both dragging by two-body forces 
and large artificial jumps in position. In Models SDH and 
BS, the artificial advection is turned off when tubu > 10m g ; 
in the fiducial (low) resolution, the motions become smaller 
(larger) and more chaotic. We note that SDH05 and BS09 



both use a multiphase star formation model, which produces 
a smoother temperature profile around the black hole than 
our star formation model. This difference may explain the 
poor advection in Models SDH and BS. 

From the simulations, the advection methods in Models 
ONB and WT appear to have some advantages; they do not 
add mass to the galaxy as in the tracer mass method, nor do 
they couple to gas particles. By associating the direction of 
motion to stellar densities (Model ONB) or the local centre 
of mass within r- in f (Model WT) , undo influence of gas parti- 
cles, which are dynamically and continuously affected by the 
AGN feedback, is avoided. Lastly, by limiting the distance 
a black hole can be artificially advected per iteration, it is 
possible for it to remain in the centre of a void. 

The advection algorithms in Models DQM, ONB and 
WT are comparatively unaffected by resolution; this is ex- 
pected since the tracer mass is independent of the local gas 
particles, and local averages vary little with resolution. The 
advection algorithm in Models SDH and BS are more reso- 
lution dependent; with lower resolution comes a larger n n f, 
thus the black hole has a larger distance it can be artifi- 
cially displaced. It is worth noting that Fig. [5] clearly shows 
that the oscillations of Model SDH is damped with higher 
resolution, but is not completely removed. 

Due to the notable black hole motion in Model BS, we 
ran a variation, Model BSw, which uses the black hole advec- 
tion algorithm of Model WT. In this model, the black hole 
had smoother accretion rates and more continuous gas den- 
sities and temperatures near the black hole. The outbursts 
near apoapsis were stronger and more symmetric between 
the two galaxies than in Model BS. Ironically, although we 
might argue that Model BSw is inherently better than Model 
BS because of the better black hole trajectories, both mod- 
els yield remnants with similar total black hole masses and 
stellar velocity dispersions, even though the black holes in 
Model BS never merge and the feedback history of the two 
models is different. 

3.2 Accretion rate 

The mass evolution of the black hole is a key aspect of the 
AGN feedback models and is an important diagnostic since 
it sets the Mbh-o" results. In the top panel of Fig. |6j we 
show the total black hole mass over time; in the bottom 
panel, we show the accretion rates (geometrically averaged 
over both black holes and plotted in bins of 10 Myr); in the 
middle panel, we show when a black hole is accreting at its 
Eddington limit. There are three major epochs of black hole 
growth: at the beginning (Models ONB, WT and DQM), at 
apoapsis (Models SDH, BS and WT), and at core merger 
(Models SDH, BS and DQM). It is typically at these times 
when a black hole is undergoing Eddington accretion; see 
Table [5] for the percent of time that a black hole is accreting 
at its Eddington limit. 

Three models use the Eddington-limited Bondi accre- 
tion rate: Models SDH, BS and WT. In SDH05 (see bottom 
panel of their fig. 14), the accretion rate increases from the 
onset of the simulation until prior to second periapsis, at 
which point there is only a slight decrease. At core merger, 
there is an increase of approximately 1.8 dex, followed by 
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Figure 6. Bottom: The accretion rates for each of the fiducial 
resolution models geometrically averaged over both black holes 
and plotted in bins of 10 Myr. Top: The total black hole mass 
in each fiducial resolution simulation. Middle: Points represent 
when a black hole is accreting at Afgdd- The vertical line of the 
same linestyle indicates when the black holes merge. The four 
black vertical lines indicate the time of first periapsis, apoapsis, 
second periapsis and core merger. 



a sudden drop as the system is totally disrupted. Since our 
initial conditions and our Model SDH was constructed to 
mimic SDH05, we expect similar results, which, within rea- 
son, are obtained. In both Models SDH and WT, there is 
an initial decline in accretion rate as the galaxies are relax- 
ing; this initial relaxation does not occur in SDH05 since 
the vertical structure of the gas is set in hydrostatic equi- 





Pre-BH Merger 


Merged BH 


Model SDH 


8.32 ( 6.82) 


1.95 (0.04) 


Model BS 


5.87 ( 2.26) 


-(-) 


Model ONB 


5.78 ( 3.17) 


0.00 (0.00) 


Model WT 


9.99 ( 5.75) 


0.00 (0.00) 


Model DQM 


2.35 ( 8.22) 


6.86 (1.42) 


Model BSw 


2.11 ( 1.18) 


0.00 (0.00) 


Model ONBc 


0.00 ( 0.37) 


0.00 (0.00) 


Model WTh 


8.78 ( 6.94) 


2.10 (0.00) 


Model DQMe 


17.72 (15.74) 


1.20 (0.00) 



Table 5. Percent of time black holes (BHs) are accreting at their 
Eddington limit; the pre-BH merger times are averaged over both 
black holes. Dashed lines indicate that the black holes did not 
merge prior to 1.5 Gyr. The percentages for fiducial resolution 
Models ONB and ONBc are taken at 1.1 Gyr and 1.25 Gyr, re- 
spectively. 



librium, thus this initial difference is expected. After first 
periapsis, gas is funnelled into the core, thus, in agreement 
with SDH05, there is an increasing accretion rate. At apoap- 
sis, there is a plateau in Model SDH, but a short-lived spike 
in Model WT. At core merger, there is an increase followed 
by a rapid decrease as the systems are totally disrupted. 

For the first 400 Myr, the accretion rate of Model BS 
follows the same trends as Models SDH and WT. During 
this time, a(nn) > 100, which explains why MModoi bs > 
Miviodois sdh & wt for 230 < i/Myr < 400. After this point, 
AfModei bs decreases, which is a result of the strong outflows 
from feedback and the motion of the black holes; the latter 
prevents the black holes from remaining in a steady environ- 
ment, resulting in rapidly varying gas characteristics around 
the black hole. The accretion rates in Model BSw continue 
to increase until apoapsis, after which they decrease due to 
strong outflows, but to a smoother accretion profile that is 
slightly higher than in Model BS; see Fig. [7] 

In DQM11, there is a large jump in M v isc followed by 
a slow decline shortly after first periapsis and again start- 
ing at second periapsis; see the bottom panel of their fig. 
1. Their initial galactic separation is approximately twice of 
ours, thus their evolution times are longer, which may ac- 
count for differences between our models and theirs. Also, 
their mass and spatial resolution are slightly better than 
in the models presented here, thus their feedback is likely 
distributed amongst fewer particles, which keeps their void 



relatively smaller (see discussion in E 3.3 1 . 

In Model DQM, there is not the accretion epoch be- 
tween first periapsis and apoapsis that occurs in DQM11; 
this is a result of both a shorter time during which gas can 
fall into the galactic cores, and the slightly larger r[ n f which 
allows for efficient feedback to maintain the void during 
non-cataclysmic events. At second periapsis, there is an in- 
crease of approximately 2.2 dex in Model DQM compared to 
DQMll's increase of approximately 3 dex. After the merger, 
both accretion rates drop off, but DQM 11 falls off faster. 

ONB08 is a cosmological simulation run from z = 49 
to z = 0. Over the final 7 Gyr of their simulation, there are 
no significant galaxy mergers and there is negligible black 
hole growth (see their fig. 8). During this time, their accre- 
tion rate lies comfortable within the RAIF regime by more 
than an order of magnitude. In Model ONBc, there is a low 
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Figure 7. The accretion rate for Models BS and BSw, geomet- 
rically averaged over both black hols and plotted in bins of 10 
Myr. The vertical line of the same linestyle indicates when the 
black holes merge. The four black vertical lines indicate the time 
of first periapsis, apoapsis, second periapsis and core merger. 



accretion rate that is continually decreasing and comfort- 
ably within the RAIF regime; however, this model uses a 
very conservative search algorithm to calculate A/drag, thus 
the low accretion rate is likely artificially low because of 
the numerical scheme. When we use a less stringent search 
algorithm that is better suited for use in Hydra at the 
resolutions we present, the accretion rate quickly rises to 
Mdrag ~ 0.005 MQyr -1 and remains approximately con- 
stant for the remainder of the simulation. This relatively 
constant rate is a result of the feedback energy being de- 
posited into the halo gas rather than modifying the galactic 
core. Since the feedback is being returned to the halo gas, 
there are negligible core differences between Models ONB 
and ONBc, except for the accretion rate. 
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Figure 8. Gas density (top) and temperature (bottom) within 
r; n f the black holes, geometrically averaged over both black holes 
and plotted in bins of 10 Myr. The vertical line of the same 
linestyle indicates when the black holes merge. 



3.3 Energy feedback algorithm 

The two broad categories of returning the feedback energy 
are thermally, E, and kinetically, p. Thermal energy is re- 
turned to gas particles in Models SDH, BS, ONB and WT by 
increasing the internal energy of the particles; kinetic energy 
is returned to gas particles in Model DQM by increasing the 
momentum of the particles. Although these two approaches 
are related by E = pc, both forms of feedback yield dras- 
tically different results, which are readily apparent on both 
large scales and small scales (e.g. see Fig. [8] for gas density 
and temperatures around the black holes). 

The feedback in Model ONB is returned to the gas par- 
ticles with p < 5 x 10 -26 g cm -3 ; thus, this energy is dis- 
tributed to halo gas particles with distances 0.5 < r/kpc < 
3.5 (7 < r/n n f < 48) from the black holes, which leaves the 
galactic cores relatively unmodified by AGN feedback. As 
gas falls towards the galactic cores (due to tidal interactions 
at first periapsis), it remains there since the AGN feedback 
is not local enough to remove it, hence these models have 
the highest core density. The major core heating events in 



this case are from shock heating of the infalling gas after 
first periapsis and from the final core collisions. 

The AGN feedback in the remaining models is deliv- 
ered directly to the core region, leading to a more distinct 
modification of the entire system. From the onset of the re- 
maining thermal feedback simulations, there is a higher core 
temperature and a lower core density than in Model ONB. 
Tidal interactions at first periapsis cause gas to be funnelled 
towards the cores, but unlike Model ONB, the temperature 
continues to rise from AGN heating, and the density remains 
moderated. Shortly after first periapsis, a galactic bar starts 
to form and persists strongly for approximately 200 Myr; 
as it dissipates there is a nuclear inflow of gas, accounting 
for the local density maxima near apoapsis in Models SDH 
and WT. After the dissipation of the bar, the galaxies be- 
gin to eject gas, preferentially along the polar axis. As the 
cores merge, there is a final set of feedback events as the gas 
is heated and blown away; after this, accretion rates and 
core densities drop while core temperature remains approx- 
imately constant. 
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Kinetic feedback affects the galactic core in different 
quantitative and qualitative ways than thermal feedback. In 
this model, the energy is returned radially to gas particles 
as a momentum boost, which quickly and efficiently forms 
a void around the black holes; in Model DQM, a void of 
r ~ 0.55n n f is formed within the first 60 Myr. Once the void 
is formed, it persists for the remainder of the simulation. 
Strong events, such as the infalling gas after first periap- 
sis or the collision at core merger, can briefly decrease the 
void radius, increasing the core density. Otherwise, the void 
and core characteristics are held constant, as we can see by 
an approximately constant temperature and density profile 
between first periapsis and core merger. 

Model DQMe uses the same kinetic feedback algorithm 
as in Model DQM, but instead r- m { = 4es2, which leads to 
an initially higher accretion rate. The increased amount of 
feedback energy is returned to more particles, delaying the 
formation of the void. When the void forms after 140 Myr, 
it has a radius of r ~ 0.85?"i n f ~ 1 kpc, which is highly unre- 
alistic. This void is more rigid than the void in Model DQM, 
only being affected at second periapsis and core merger. The 
unrealistically large r- m f in Model DQMe also leads to a pre- 
mature merger of the black holes, which is poorly handled 
by the merger subroutine. Thus, the importance of choosing 
a physical r in f is clear: catastrophic and unphysical results 
will ensue if it is not chosen appropriately. 

3.4 Particle accretion algorithm 

A particle accretion algorithm should accrete particles such 
that the internal mass, Mbh, and dynamical mass, tubh, of 
any given black hole remain similar for all time. In Model 
DQM, a tracer mass is used to represent the black hole parti- 
cle, rendering the comparison between internal and dynamic 
masses moot; these models use a continual-conditional al- 
gorithm to remove 'accreted' gas, thus any conclusions we 
make regarding this category can be applied. In the re- 
mainder of our fiducial resolution models, the first accretion 
event will increase the dynamical mass by 36 per cent since 
fttBH = 0.36m g . All subsequent accretion events will add 
relatively less mass to the black hole particle. Thus we ex- 
pect the ratio Am /Mbh, where Am = wbe — Mbh, to start 
at |Am/MeH| < 0.36 and decay to zero with time. In the 
top panel of Fig. |9j we plot Am/ Mbh for one black hole in 
Models SDH, BS, ONB and WT. The expected behaviour is 
obtained in Models BS, ONB and WT; the decline in Mod- 
els ONB and WT is smoother than in Model BS since there 
is no stochastic component to the accretion. The expected 
decline is not observed in Model SDH, where the ratio has 
no apparent trend, and at late times Am/A/sH — > —0.01 
rather than zero. 

The black holes do not all grow at the same rate. Thus, 
for a normalised comparison, we plot Am/m g in the bottom 
four panels of Fig. [9] for Models SDH, BS, ONB and WT, 
respectively. Since every accretion event adds mass m g or 
m.g/2 to the dynamical mass, we expect |Am/m g | 1 for 
all time, unless multiple particles are accreted in one step, 
which can occur in the stochastic models. 

In Model SDH, particle accretion events are only depen- 
dent on the environment. This lack of Mbh-tibh coupling 
permits the internal and dynamic masses to follow different 
histories, as seen by the lack of trend in the second panel 
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Figure 9. A comparison of internal mass, Mbh, an d dynamical 
mass, rriBH, for one black hole in selected fiducial runs. Top: The 
ratio Atti/Mbhi where Am = tiibh — Mbh- Bottom Four Panels 
(top to bottom): The ratio Am/m g for Model SDH (stochastic- 
unconditional), Model BS (stochastic-conditional), Model ONB 
(continual-conditional) and Model WT (continual-conditional 
particle accretion). The vertical line of the same linestyle in the 
top panel indicates when the black holes merge. 



of Fig. [9] The discontinuity in Am/m g at core merger is a 
result of both black holes following different mass growth 
histories, both internally and dynamically. 

The stochastic-conditional accretion algorithm in 
Model BS allows for fluctuations in an otherwise smooth 
process, and the probability's dependence on Am forces 
good agreement between the internal and dynamic masses. 
Finally, by design, Model ONB (WT) always maintains 
0.0 < Am/m g < 1.0 (-0.5 < Am/m g < 0.5), displaying 
the effectiveness of the continual-conditional particle accre- 
tion methods. 
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Table 6. Merger time of the black holes, as measured from the 
beginning of the simulation. Dashed lines indicate that the black 
holes did not merged within our 1.5 Gyr simulation time. 

3.5 Black hole merger algorithm 

All of the models include a black hole merger prescription, 
which conserves internal, dynamical and tracer mass; see 
Table [6] for precise merger times. Thus, the merger should 
have limited gravitational impact on the surrounding envi- 
ronment. However, the amount of feedback energy returned 
pre- and post-merger will vary depending on the model, 
which we analyse below. 



Consider two identical black holes in the merged core 
during the black hole merger. First, the total region of influ- 
ence of the unmerged black holes will be between V = f 7rrf nf 
and 2V depending on separation of the black holes. Since 
n n f is independent of black hole mass, the total region of 
influence after the merger will simply be V . In all cases, 
E oc M, so the amount of feedback energy available prior 
to the merger will be proportional to Mbh ± + Mbh 2 after 
the merger, we will have energy proportional to Mbh 1 +bh 2 ', 
if we define M to be the accretion rate of one black hole 
prior to the merger, then for each accretion rate we studied, 
across the merger we have 

2A/Edd — > 2A/Edd, 

2Mb -> 4A4b , 

2Mdrag — > Afdrag, 
2M viac -> M visc . 

We can clearly see that the only total accretion rate in the 
core that remains unchanged by the merger is A^Edd! the 
total Bondi accretion rate doubles and the total viscous and 
drag accretion rates halves. Thus, in Models SDH, BS and 
WT (which use Bondi accretion), there is more feedback 
energy available after the merger than before, and this in- 
creased amount of energy will be distributed to a smaller 
region. Likewise, there is less feedback energy available in 
Models ONB and DQM after the merger, which will also 
be distributed within a smaller volume. Thus, we can ar- 
gue that the galactic core becomes more energetic assuming 
Bondi accretion, but less energetic assuming viscous or drag 
accretion when two black holes merge. 

In the models, Model DQM has the least stringent 
merger prescription, with only a separation criteria. As ex- 
pected, this is the first model to merge. The most stringent 
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Figure 10. Global star formation rate for the fiducial (top) and 
low resolution (bottom) runs. The vertical lines represent the 
times of first periapsis, apoapsis, second periapsis and core merger 
for the dynamic black hole mass models, as given in Table [4] 



merging criteria is in Model BS, which is never met. The 
circular velocity at the more massive black hole's smoothing 
length is typically ten times smaller than the local sound 
speed, and the relative velocity between the black holes 
never drops below this threshold. Moreover, the chaotic mo- 
tion and two-body interactions of the black holes typically 
keeps them further apart than the separation criteria. 



3.6 Global star formation rates 

During a major merger, there is typically an burst of star for- 
mation at apoapsis and again during core merger due to the 
infall of gas on to the galactic core (e.g. |Mihos fc Hernquist 
(19961); SDH05; DQM11). There are several characteristics 



that can modify the star formation rate at these epochs, 
including initial conditions and the inclusion of a galactic 
bulge and/or AGN feedback (SDH05). In Hydra, the star 
formation algorithm differs from those used in SDH05, BS09, 
ONB08 and DQM11, so we do not expect an exact repro- 
duction of their results; however, any differences in SFRs 
amongst our models will be solely the result of the AGN 
feedback algorithm and its interaction with the star forma- 
tion algorithm. The time-averaged SFR for the fiducial (top 
panel) and the low resolution (bottom panel) are plotted in 
Fig. [TO] 

In these models, there are notable star formation bursts 
shortly after the beginning of the simulation, at apoapsis 
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and again at core merger. In all cases, the initial burst is 



a result of the gas disc relaxing (Williamson & Thacker 
|2012h; this initial burst also exists in simulations of an iso- 



lated galaxy that is constructed in a similar manner to the 
one presented here. SDH05 and DQM11 create the vertical 
structure of their gas disc by assuming hydrostatic equilib- 
rium after setting their equation of state. Although cooling 
is not considered during this process, it means that their 
initial conditions and simulations are more equivalent than 
ours. Thus, as expected, this initial burst of star formation 
is missing from their results. This initial burst is resolution 
dependent, being higher and shorter lived in the fiducial res- 
olution models. After this epoch in the low resolution mod- 
els, the only outstanding features are post-apoapsis burst 
in Model SDHi, the core merger burst for Models SDHi and 
DQMi, and the expected drop in SFRs for all models (except 
Model ONBi) shortly after core merger. The low resolution 
features do not necessarily have counterparts in the fiducial 
resolution simulations and vice versa, as a result of the stel- 
lar feedback being resolution dependent: low resolution star 
formation events extend over a considerably larger volume 
and mass than fiducial resolution events due to the larger 
gas particle mass and lower spatial resolution. Thus, a low 
resolution event has an unavoidably greater impact on its 
environment, and causes comparatively large regions of the 
galaxy to undergo a cooling delay. This explains the delayed 
star formation peak after apoapsis in Model SDHi, the lack 
of peak in Model WTi after apoapsis, and generally lower 
SFRs at core merger. 

In Model SDH, the SFR increases by 0.9 dex at both 
apoapsis and core merger; in SDH05, there is no burst at 
apoapsis and there is an increase of 1 dex at core merger. 
This burst at core merger exists in all the models, except 
Model WT, where the high feedback rate of thermal energy 
prevents the gas from cooling in to the regime where star 
formation can commence. After the core merger burst, the 
SFRs drop, as in SDH05, due to the heating and total dis- 
ruption of the system. 

In Model WTh, the star formation burst at apoapsis 
peaks approximately 50 Myr earlier than in Model WT. 
The apoapsis burst is a result of a nuclear inflow of gas, 
and the lower value of /i m i n allows higher densities to be re- 
solved sooner; thus inflowing gas enters the permissible den- 
sity regime of star formation for Model WTh before Model 
WT. For the same resolution argument, Model WTh also 
has a small burst of star formation (an increase of 0.6 dex) 
at core merger. This difference in the SFR highlights one of 
the well known drawbacks of resolution dependency in the 
star formation algorithm, although developing a model that 
is independent of resolution awaits a full understanding of 
star formation in a global context. 

The SFR in Model ONB remains relatively smooth until 
a small increase of 0.6 dex at core merger. Even after core 
merger, the SFR remains relatively unchanged due high core 
densities and moderate core temperatures. 

In DQM11, there is a major star formation burst shortly 
after first periapsis and again starting at second periapsis 
(see their fig. 1, middle panel); Model DQM reproduced the 
burst at second periapsis. The hot ring around the void 
means that the star formation must occur in the galactic 
disc, which results in the lack of peak after first periapsis. 
At second periapsis, the merger is efficient enough to cre- 



ate a galactic core, causing the core merger star formation 
burst. After this final burst, the SFR drops to small values, 
agreeing with DQM11. 



4 FINAL STATES 

We evolved our galaxy merger simulations for 1.5 Gyr, and 
the characteristics of the remnants are presented here. The 
exception is the fiducial resolution Model ONB, which was 
evolved for 1.1 Gyr. From our analysis of Model ONBi, we do 
not expect any significant evolution between 1.1 and 1.5 Gyr. 
Therefore, for the final remnant of Model ONB, we present 
a combination of data from 1.1 Gyr, data extrapolated to 
1.5 Gyr, and data from Model ONBi; we will clearly state 
where the data comes from in each instance. 

The final remnant in all models is a reformed gas and 
stellar disc surrounded by a fragmented hot gas halo. The 
radial profiles for gas temperature, gas density and stellar 
column density (averaged over all azimuthal and polar an- 
gles) for each remnant are similar, as seen in Fig. 
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the 

values typically span 0.5 dex at any given radius. In Figs. 
|12| and |13| we have plotted the final gas temperature, gas 
column density and stellar column density for frames of 100 
and 20 kpc per side, respectively. 



4.1 Gas properties 

By the end of the simulation, the gas that was originally in 
a galactic disc is heavily depleted while the majority of the 
gas that was originally in the halo remains; the fraction of 
remaining gas for each fiducial resolution model is given in 
Table[7| The majority of the depleted gas has been converted 
in to stars as opposed to being accreted on to a black hole. 
The gas that is ejected from the disc has temperatures and 
pressures a few orders of magnitude greater than the halo 
gas, thus it can easily escape from the system. Therefore, 
the gas halo does not play a roll in recycling the ejected gas 
in to stars. 



The amount of substructure varies considerably 
amongst remnants. Model ONB's remnant is a well formed 
disc in a smooth halo, while Models BS and WT have a 
loosely reformed disc in a fragmented halo. Shortly after 
core merger, Model BS undergoes an outburst event which 
expels most of the gas, and by 1.5 Gyr, this gas is starting to 
recondense on to the core. In Model WT, the remaining gas 
cools and fragments, and then begins to recondense on to 
the core. Thus, in both models there are similar remnants at 
1.5 Gyr but different histories prior to this time. In all cases, 
a small disc begins to reform; the surface density profile for 
the reformed fiducial resolution discs is plotted in Fig. |14| 
The scale lengths of the discs embedded within the stellar 
remnant are approximately 0.5 kpc, compared to the initial 
value of 2.46 kpc. There are voids in the centre of the discs 
in Models SDH, BS and WT, thus the peak surface density 
occurs at 0.7-1.5 kpc from the centre. The voids persist from 
the residual angular moment of the gas and from the weak 
(but remaining) AGN feedback. 




Figure 12. Left to right: Face-on gas temperature, gas column density, and stellar column density of each remnant. Each frame measures 
100 kpc per side, with a image resolution of 195 pc pixel -1 (391 pc pixel -1 ) for the fiducial (low) resolution models. 
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Figure 13. Left to right: Face-on and edge-on gas temperature, face-on and edge-on gas column density, and face-on and edge-on stellar 
column density of each remnant. Each frame measures 20 kpc per side, with a image resolution of 39 pc pixel -1 (78 pc pixel -1 ) for the 
fiducial (low) resolution models. 



4.2 Stellar remnant 

For each remnant, we calculate the stellar velocity disper- 
sion, a, around the black hole. To mimic observations, we 
calculate this using 



Jo °<r loB (r)I(r)rdr 
f* e I(r)rdr 



(25) 



where a\ os is the line of sight velocity dispersion, I(r) is the 
projected 2D stellar density profile, and R e is the half-light 
(mass) radius; this is similar to what is done numerically 



in DQM11 and observationally in Giiltekin et al. (20091 



Since the final stellar configuration is highly triaxial (see 
the two right-most columns in Fig. \13\ , the a we present in 
Fig. [15] is averaged over fOOO random lines of sight; we have 



also included the full range of <r's calculated, as well as the 
preferentially chosen lines of sight along the ±x— , ±y— and 
±z— directions. 

For the fiducial (low) resolution models, the velocity 
dispersion spans a range of ~75 ± 15 km s -1 (~90 ± 15 km 
s _1 ), and is latitude dependent. Since the remnant reforms 
a disc, there is a large a if the line of sight lies in the plane of 
the disc, and a decreases as the elevation increases. Except 
for Models BS and ONB, the average a lies is within the 
one-sigma scatter of the observational relation. 

In Model BS, the black holes never merge and do not 
fall within one-sigma of the observed relation; the black holes 
in Model BSw merge, but still do not fall within one-sigma 
of the observed relation. This indicates that the difference 
from observational expectations is not a result of the lack of 
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0.988 
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ONB 


1.053 
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40.6 


0.00165 


WT 


1.061 
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0.452 
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0.00450 


DQM 


1.098 


0.577 


0.127 


0.987 


165.9 


0.00339 


BSw 


1.046 


0.802 


0.592 


0.992 


16.6 


0.00069 


ONBc 


1.055 


0.766 


0.513 


0.996 


1.4 


0.00000 


WTh 


1.065 


0.716 


0.415 


0.991 


169.8 


0.00521 


DQMc 


1.064 


0.720 


0.423 


0.991 


615.1 


0.01917 



Table 7. Ratios of final to initial global mass components for each of our fiducial resolution models; black hole masses are dynamic 
masses (except for Model DQM) since gas is removed in discrete amounts of m gas . Gas labelled as disc gas or halo gas is based upon its 
initial position. We define M acc as the total mass of gas (dynamically) accreted by black holes and AAfg aa = Af gas ; — M gas f. Depleted 
gas that is not accreted on to a black hole has been converted into stars. The values for Models ONB and ONBc are extrapolated to 
1.5 Gyr based upon data from the final output and their low resolution counterparts. Specifically, we assume an average SFR of 1.0 
Mq yr _1 (1.6 Mq yr _1 ) and M B H = 5.6 X 10~ 3 M Q yr _1 (6.0 X 10~ 8 Mq yr _1 ) for the remainder of the simulation for Model ONB 
(ONBc). 
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Figure 15. The Mbh-t relation for our simulations, along with the observed relation (red) and the one-sigma scatter (dashed) from 
|Giiltekin et al.| ( |2009| ). For our five primary models at each resolution, the dot represents the average a of 1000 random lines of sight, 
the horizontal bars represent the range of all calculated it's, and the remaining three symbols on the horizontal line represent a taken 
along the ±x— , iy— and ±2— lines of sight. For our variant models, we have only plotted the average a of 1000 random lines of sight 
(solid triangles). We have explicitly labelled which points/set of points corresponds to which models. The black dots in the red rectangles 
represent the black holes in Model BS; the lower two points represent the actual black holes, and the upper point is calculated assuming 
one black hole with Mbh = ^BHx + Mbh 2 at the centre of mass of the black hole system. The solid triangle in the red rectangle is the 
result for Model BSw. The points in the green box are our initial relationships. The data for Models ONB and ONBc are taken at 1.1 
and 1.25 Gyr, respectively. 



merger for the model parameters we have implemented. In 
these models, between apoapsis and second periapsis, the ac- 
cretion rates are lower than the other models, which hinders 
black hole growth. Also, randomly depositing feedback en- 
ergy around the black hole, rather than returning it isotrop- 
ically or kernel- weighted, could increase the velocity disper- 
sion of the stars by increasing the anisotropy of the gas. 

In Model ONB, the black hole mass is 9.5 times lower 
than predicted by the Mbh-o" relation, assuming the stellar 
velocity dispersion remains constant; a mass 3.4 times larger 



would bring the point within the one-sigma scatter of the 
relation. Since the feedback energy is returned far from the 
galactic core, gravity is the only influence the black hole 
has on the nearby stars. Thus, the self-regulating feedback 
mechanism is never fully realised, which could explain the 
lack of agreement with the observed relation. This is also 
true for Model ONBc, but to a more extreme extent since the 
low accretion rate means that the black hole only negligibly 
grows and cannot couple to the stellar system. 
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Figure 11. Radial profiles of our remnants, averaged over all az- 
imuthal and polar angles. From top to bottom, we show gas tem- 
perature, gas density, and stellar density. The profile for Model 
ONB is taken at 1.1 Gyr. 

5 CONCLUSIONS 

This paper details a study of five different AGN feedback al- 
gorithms, which were all run using the numerical code Hy- 
dra, started from the same initial conditions, and imple- 
mented the same star formation algorithm; two resolutions 
were tested for each model. The AGN feedback algorithms 
used in Models SDH, BS, ONB and DQM were previously 
found in the literature and the algorithm used in Model WT 
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Figure 14. Gas surface density profile, averaged over all az- 
imuthal angles, for the initial and final discs of the fiducial reso- 
lution models. The profiles are truncated at the edge of the disc, 
and plotted in bins of 49 pc. The remnant profile for Model ONB 
is taken at 1.1 Gyr. 

was developed specifically for this study. We also ran four 
additional models, each of which was a slight variation of 
one of the primary models to study a specific aspect of the 
model. 

We have performed a full analysis of all models, but 
only presented fiducial resolution results of the five primary 
models and selected other results. Although the low resolu- 
tion models run considerably faster, their results are noisy, 
and some of the physics is essentially lost. However, the low 
resolution studies are a good test of the varying feedback 
models, and provide a first-order comparison between the 
models. Furthermore, the continued implementation of cos- 
mological simulations mean these resolutions will inevitably 
be used in the future. 

By comparing these models using the same initial con- 
ditions, star formation algorithm and numerical code, any 
differences obtained are a direct result of the AGN feed- 
back. Although we have tried to isolate each component of 
the AGN feedback algorithm, we accept that they are all 
intimately intertwined, and the effects cannot necessarily 
be disentangled from one another. We also accept that ev- 
ery simulation we run involves free parameters; although we 
have - as best as possible - matched them to their source 
simulations from the literature, it is possible that we can 
modify these parameters such that all simulations will re- 
turn similar final states. Since the goal of this paper was to 
compare different feedback models and not to analyse any 
given model in great detail, we defer to the source litera- 
ture for discussions on the sensitivity of the results to the 
free parameters. We thus refrain from making strong state- 
ments about the 'accuracy' of given models. Our principal 
conclusions of the five key components are as follows: 

(i) The black hole advection algorithm plays a key role 
since small displacements can cause great changes in the ac- 
cretion rate. Coupling the black hole to a gas particle can 
allow for oscillations or chaotic motion if a void is formed 
around the black hole, or if the nearby gas does not meet the 
velocity criteria. Using a tracer mass yields smooth move- 
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ment of the black hole, but the evolution time of the merger 
is decreased and disc morphology is altered in non-trivial 
ways such as preventing the formation of a bar. This obvi- 
ously impacts other physical properties, most notably the 
SFR as well. Clearly evolving the mass of the black hole in 
a way that does not overly impact the expected evolution 
of the disc is an important goal. The algorithms that limit 
the distance a black hole can be displaced per iteration and 
where the direction is based upon the local stellar or total 
mass distribution appear to be optimal. 

(ii) We test four different accretion rates: Mb (a = 100), 
Mb (a = a(nn)), Md rag and Af v i sc . Each model yields a dif- 
ferent accretion history, with substantially different qualita- 
tive profiles and quantitative differences of up to a factor 
of three orders of magnitude at any given time. The total 
black hole mass in the remnants varies by factor of 6.9, with 
Model BS being the least massive, and DQM being the most 
massive. The three models that fall within one-sigma scat- 
ter of the Mbh-o" relation have final black hole masses that 
differ by less than 30 per cent. 

(iii) The form of feedback can drastically modify the 
large- and small-scale systems. Delivering the feedback to 
the halo gas leads to little modification of the core region, 
resulting in high core densities and nominal temperatures. 
Thermal feedback delivered to the core region can drive 
strong outflows. Kinetic feedback delivered to the core region 
efficiently creates a void around the black hole, resulting in 
low gas temperature and density. This feedback is persistent 
and efficient, keeping the core properties nearly constant for 
all time; only a cataclysmic event, such as a core merger, 
can modify these characteristics. 

(iv) We tested three categories of particle accretion algo- 
rithms: stochastic-unconditional, stochastic-conditional and 
continual-conditional; two of the three continual-conditional 
algorithms randomly selected which particles were accreted 
but never how many or how often. The continual-conditional 
cases gives the best agreement between the dynamical and 
internal mass; here, the discrepancy is never more than 2m g . 
The stochastic- unconditional algorithm of Model SDH, how- 
ever, contains discrepancies up to 15m g . Thus, for agree- 
ment between dynamical and internal masses, continual- 
conditional algorithms appear optimal. 

(v) In Models BS and BSi, the black holes never merge; in 
these models, there is considerable chaotic movement of the 
black holes in the remnant core, preventing both merger cri- 
teria from being simultaneously met. In the remainder of the 
models, the black holes merge during or shortly after core 
merger, as one would reasonably expect. The importance 
of a 'reasonable' merger time is that the amount of feed- 
back energy available increases across the merger for Bondi 
accretion and decreases for drag and viscous accretion. 

The models presented here do not represent an exhaus- 
tive list. For example, an accretion rate suggested by |Hobbs] 
|et al.| ( |2012| ) is an interpolation between Bondi accretion and 
a free-fall accretion rate. Power et al. (20111 suggest using 



an accretion disc particle; in this two-tier accretion process, 
gas is first accreted on to an accretion disc, and then the gas 
is accreted from the disc on to the black hole using a viscous 
time-scale. The accretion disc particle method of accretion 
has recently been implemented in major merger simulations 
similar to those presented here ( Wurster fc Thacker||2013 1. 



Lastly, while it is possible that free parameters in differ- 
ent models can be adjusted such that all yield similar final 
states, it is very clear that this will not produce agreement 
throughout the evolution. While we have a vast amount of 
observational data of many stages of major mergers, includ- 
ing final remnants, it is difficult to build these various snap- 
shots into an obvious picture of the evolutionary processes, 
and the fundamental theoretical development of AGN mod- 
elling still has much to contribute in this regard. As we have 
shown here, given the many different models of AGN feed- 
back, building an understanding of the evolution of mergers 
is an important and challenging task. 
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